Wave - equation migration velocity analysis . I . Theory

We present a migration velocity analysis (MVA) method based on wavefield extrapolation. Similarly to conventional MVA, our method aims at iteratively improving the quality of the migrated image, as measured by the flatness of angle-domain commonimage gathers (ADCIGs) over the aperture-angle axis. However, instead of inverting the depth errors measured in ADCIGs using ray-based tomography, we invert ‘image perturbations’ using a linearized wave-equation operator. This operator relates perturbations of the migrated image to perturbations of the migration velocity. We use prestack Stolt residual migration to define the image perturbations that maximize the focusing and flatness of ADCIGs. Our linearized operator relates slowness perturbations to image perturbations, based on a truncation of the Born scattering series to the first-order term. To avoid divergence of the inversion procedure when the velocity perturbations are too large for Born linearization of the wave equation, we do not invert directly the image perturbations obtained by residual migration, but a linearized version of the image perturbations. The linearized image perturbations are computed by a linearized prestack residual migration operator applied to the background image. We use numerical examples to illustrate how the backprojection of the linearized image perturbations, i.e. the gradient of our objective function, is well behaved, even in cases when backprojection of the original image perturbations would mislead the inversion and take it in the wrong direction. We demonstrate with simple synthetic examples that our method converges even when the initial velocity model is far from correct. In a companion paper, we illustrate the full potential of our method for estimating velocity anomalies under complex salt bodies. I N T R O D U C T I O N Seismic imaging is a two-step process: velocity estimation and migration. As the velocity function becomes more complex, the two steps become increasingly interdependent. In complex depth imaging problems, velocity estimation and migration are applied iteratively in a loop. To ensure that this iterative imaging process converges to a satisfactory model, it is crucial that the migration and the velocity estimation are consistent with each other. Paper presented at the EAGE/SEG Summer Research Workshop, Trieste, Italy, August/September 2003. ∗E-mail: [email protected] Kirchhoff migration often fails in areas of complex geology, such as subsalt, because the wavefield is severely distorted by lateral velocity variations leading to complex multipathing. As the shortcomings of Kirchhoff migration have become apparent (O’Brien and Etgen 1998), there has been renewed interest in wave-equation migration, and computationally efficient 3D prestack depth migration methods have been developed (Biondi and Palacharla 1996; Biondi 1997; Mosher, Foster and Hassanzadeh 1997). However, no corresponding progress has been made in the development of migration velocity analysis (MVA) methods based on the wave equation. We aim at filling this gap by presenting a method that, at least in principle, can be used in conjunction with any downward-continuation migration method. In particular, we have been applying our new C © 2004 European Association of Geoscientists & Engineers 593 594 P. Sava and B. Biondi methodology to downward continuation based on the double square root (Yilmaz 1979; Claerbout 1985; Popovici 1996) or common-azimuth (Biondi and Palacharla 1996) equations. As in the case of migration, wave-equation MVA (WEMVA) is intrinsically more robust than ray-based MVA because it avoids the well-known instability problems that rays encounter when the velocity model is complex and has sharp boundaries. The transmission component of finite-frequency wave propagation is mostly sensitive to the smooth variations in the velocity model. Consequently, WEMVA produces smooth, stable velocity updates. In most cases, no smoothing constraints are needed to ensure stability in the inversion. In contrast, ray-based methods require strong smoothing constraints to avoid divergence. These smoothing constraints often reduce the resolution of the inversion that would otherwise be possible, given the characteristics of the data (e.g. geometry, frequency content, signal-to-noise ratio, etc.). Eliminating, or substantially reducing, the amount of smoothing increases the resolution of the final velocity model. A well-known limitation of wave-equation tomography or MVA is represented by the linearization of the wave equation based on truncation of the Born scattering series to the first-order term. This linearization is hereafter referred to as the Born approximation. If the phase differences between the modelled and recorded wavefields are larger than a fraction of the wavelet, then the assumptions made under the Born approximation are violated and the velocity inversion methods diverge (Woodward 1992; Pratt 1999; Dahlen, Hung and Nolet 2000; Hung, Dahlen and Nolet 2000). Overcoming these limitations is crucial for a practical MVA tool. This goal is easier to accomplish with methods that optimize an objective function that is defined in the image space (e.g. differential semblance optimization and our WEMVA) than with methods that optimize an objective function that is defined in the data space. Our method employs the Born approximation to linearize the relationship between the velocity model and the image. However, we ‘manipulate’ the image perturbations to ensure that they are consistent with the Born approximation, and we replace the image perturbations with their linearized counterparts. We compute image perturbations by analytically linearizing our image-enhancement operator (e.g. prestack residual migration) and applying this linearized operator to the background image. Therefore, the linearized image perturbations are approximations to the non-linear image perturbations that are caused by arbitrary changes of the velocity model. Since we linearize both operators (migration and residual migration) with respect to the amplitude of the images, the resulting linear operators are consistent with each other. Therefore, the inverse problem converges for a wider range of velocity anomalies than the one implied by the Born approximation. Our method is more similar to conventional MVA than other proposed wave-equation methods for estimating the background velocity model (Noble, Lindgren and Tarantola 1991; Bunks et al. 1995; Forgues, Scala and Pratt 1998) because it maximizes the migrated image quality instead of matching the recorded data directly. We define the quality of the migrated image by the flatness of the migrated angledomain common -image gathers (ADCIGs) along the apertureangle axis (Sava and Fomel 2003). In this respect, our method is related to differential semblance optimization (DSO) (Symes and Carazzone 1991; Shen 2003) and multiple migration fitting (Chavent and Jacewitz 1995). With respect to DSO, our method has the advantage that at each iteration it optimizes an objective function that rewards flatness in the ADCIGs globally (for all the angles at the same time), and not just locally as DSO does (minimizing the discrepancies between the image at each angle and the image at the adjacent angles). We suggest that this characteristic should speed up the convergence, although we have no formal proof of our assertion. This paper describes the theoretical foundations of wave-equation MVA with simple examples illustrating the main concepts and techniques. In a companion paper (Sava and Biondi 2004), we present an application of wave-equation MVA to the challenging problem of velocity estimation under salt. Here, we begin by discussing wavefield scattering in the context of one-way wavefield extrapolation methods. Next, we introduce the objective function for optimization and finally, we address the limitations introduced by the Born approximation. Two appendices detail the wave-equation MVA process and the computation of linearized image perturbations. R E C U R S I V E WAV E F I E L D E X T R A P O L AT I O N Imaging by wavefield extrapolation (WE) is based on recursive continuation of the wavefields U from a given depth level to the next by means of an extrapolation operator E, i.e. Uz+ z = Ez [Uz] . (1) Here and hereafter, we use the following notation conventions: A[x] denotes operator A applied to x, and f (x) denotes function f of argument x. The subscripts z and z + z indicate quantities corresponding to the depth levels z and z + z, respectively. C © 2004 European Association of Geoscientists & Engineers, Geophysical Prospecting, 52, 593–606 Wave-equation migration velocity analysis I 595 The recursive equation (1) can also be explicitly written in matrix form as  1 0 0 · · · 0 0 −E0 1 0 · · · 0 0 0 −E1 1 · · · 0 0 .. .. .. .. .. .. 0 0 0 · · · −En−1 1   U0 U1 U2 .. Un  =  D0 0 0 .. 0  or in a more compact notation as